library(tidyverse)
library(ggplot2)
library(haven)
library(sjPlot)
library(stargazer)
library(ggdist)
library(cregg)
library(patchwork)
library(viridisLite)
library(RColorBrewer)
colors <- viridis(option = "magma", 
                  begin = 0.1, 
                  end = .8, 
                  direction = -1, 
                  n = 9)
colorsgay<- c("#CD66FF", "#FF6599", "#FF0000", "#FF8E00", "#FFCC00", "#008E00", "#00C0C0", "#400098", "#8E008E")
df <- read_dta("data/conjoint_analysisfile.dta")


df <-df%>% 
  mutate(finance= as.factor(finance),
         elig= as.factor(elig),
         funding= as.factor(funding),
         resident= as.factor(resident),
         profile= as.factor(profile),
         legal= as.factor(legal),
         parents= as.factor(parents),
         party= as.factor(party),
         endorse= as.factor(endorse),
         IDvar= as.factor(IDvar),
         queer= as.factor(queer),
         sex= as.factor(sex),
         noncis= as.factor(noncis),
         LGB= as.factor(LGB),
         ethnic= as.factor(ethnic),
         votemain= as.factor(votemain),
         votetory= as.factor(votetory),
         ideo= as.factor(ideo),
         child= as.factor(child),
         task= as.factor(task),
         degree= as.factor(degree))

colorX<-c("#205C8A", "#E3256B")

conjoint1<- lm(selected ~ finance + elig + funding + resident + profile + legal + parents + party + endorse, data=df)
summary(conjoint1)


df$finance <- factor(df$finance,
                     levels = c(0,1,2,3),
                     labels = c("Zero compensation", "Expenses only (£12K-£18K)", "Expenses + moderate compensation (£18K-£40K)", "Expenses + large compensation (£40K-£60K)"))
df$elig <- factor(df$elig,
                  levels = c(0,1,2,3,4),
                  labels = c("Married hetero. couples", "Married couples (including homo.)", "Non-married couples (including homo.)", "Non-married hetero. couples", "Single individuals"))
df$funding <- factor(df$funding,
                     levels = c(0,1,2),
                     labels = c("No public funding", "Full funded by NHS", "Partially funded by NHS"))
df$resident <- factor(df$resident,
                      levels = c(0,1,2),
                      labels = c("Available to legal residents only", "Available to British nationals only", "Available to non-residents"))
df$profile <- factor(df$profile,
                     levels = c(0,1,2),
                     labels = c("Must be legal resident", "Must be British national", "Can reside overseas"))
df$legal <- factor(df$legal,
                   levels = c(0,1,2),
                   labels = c("Non-binding agreement", "Binding legal contract", "No legal agreement"))
df$parents <- factor(df$parents,
                     levels = c(0,1,2,3),
                     labels = c("after 3 months", "after 1 month", "after 6 months", "from birth"))
df$party <- factor(df$party,
                   levels = c(0,1,2,3,4),
                   labels = c("Conservatives", "Cross-party support", "Labour", "Lib Dems", "Greens"))
df$endorse <- factor(df$endorse,
                     levels = c(0,1,2,3,4),
                     labels = c("None", "British Fertility Society", "Stonewall", "the Economist", "Women's Aid"))
df$queer <- factor(df$queer,
                   levels = c(0,1),
                   labels = c("Hetero", "LGBT+"))
df$sex <- factor(df$sex,
                 levels = c(0,1),
                 labels = c("Man", "Woman"))
df$degree <- factor(df$degree,
                    levels = c(0,1),
                    labels = c("No degree", "Degree"))
df$noncis <- factor(df$noncis,
                    levels = c(0,1),
                    labels = c("Cisgender", "Not cisgender"))
df$votemain <- factor(df$votemain,
                      levels = c(0,1,2,3),
                      labels = c("Labour", "Conservative", "Other", "Abstain"))
df$votetory <- factor(df$votetory,
                      levels = c(0,1),
                      labels = c("Labour", "Conservative"))
df$ideo <- factor(df$ideo,
                  levels = c(0,1,2),
                  labels = c("Left", "Centre", "Right"))

df2<-df
df2 <- subset(df, child!=2) 
df2$child <- factor(df2$child,
                    levels = c(0,1),
                    labels = c("Not parent", "Parent"))



labels <- c(finance = "Financial\ncompensation", elig = "Parental\neligibility",
            funding = "Public\nfunding", resident = "Parental\nresidency",
            profile = "Surrogate\nresidency", legal = "Legal\nenforcement",
            parents = "Time of\nparentship", party = "Proposing\nparty",
            endorse = "Endorsements")

###FIGURE 4###
amce(df, selected ~ finance + elig + funding + resident + profile + 
       legal + parents + party + endorse, id = ~ IDvar)

amce.1 <- amce(df, selected ~ finance + elig + funding + resident + profile + 
                 legal + parents + party + endorse, id = ~ IDvar)


amceplot<- plot(amce.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L,
                      scales = "free_y", strip.position = "right") + 
  theme_ggdist()+
  scale_color_manual(values = colorsgay)+
  theme(legend.position = "none") +
  geom_point(size=1.5)+
  theme(strip.text = element_blank(),
        axis.title.x=element_text(face="bold"),
        axis.text.y = element_text(face="bold"))

mm.1 <- mm(df, selected ~ finance + elig + funding + resident + profile + 
             legal + parents + party + endorse, id = ~ IDvar, h0 = 0.5)

mmplot<- plot(mm.1, feature_headers = FALSE) +
  ggplot2::facet_wrap(~feature, ncol = 1L, labeller=labeller(feature = labels),
                      scales = "free_y", strip.position = "right") +
  theme_ggdist()+
  scale_color_manual(values = colorsgay)+
  geom_label(aes(label=sprintf("%0.2f", round(estimate, digits = 2))), nudge_x = -.05, size=2.5)+
  theme(legend.position = "none") +
  geom_point(size=1.5)+
  xlim(.3,.6)+
  theme(axis.title.y =element_blank(),
        strip.text = element_text(face="bold"),
        legend.position ="none",
        legend.text = element_text(face="bold"),
        plot.subtitle = element_text(face="bold"),
        axis.title.x=element_text(face="bold"),
        legend.background = element_blank(),
        legend.key = element_blank(),
        legend.title=element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.y = element_blank()) +
  geom_vline(xintercept = 0.50, linetype="dashed", colour="gray33")
amceplot+mmplot
ggsave("figures/Figure4.eps", height=24, width=20, units="cm", dpi=900)

###FIGURE 5###
mm.diff.LGB <- mm_diffs(df, selected ~ finance + elig + funding + resident + profile + 
                          legal + parents + party + endorse, 
                        ~ queer, id = ~ IDvar)

mm.diff.gender <- mm_diffs(df, selected ~ finance + elig + funding + resident + profile + 
                             legal + parents + party + endorse, 
                           ~ sex, id = ~ IDvar)

mm.diff.child<- mm_diffs(df2, selected ~ finance + elig + funding + resident + profile + 
                           legal + parents + party + endorse, 
                         ~ child, id = ~ IDvar)

mm.diff.party<- mm_diffs(df, selected ~ finance + elig + funding + resident + profile + 
                           legal + parents + party + endorse, 
                         ~ votetory, id = ~ IDvar)

mm.diff.gender<- mm.diff.gender%>% 
  mutate(statsig= if_else(p < 0.05, "Yes", "No"))

mm.diff.LGB<- mm.diff.LGB%>% 
  mutate(statsig= if_else(p < 0.05, "Yes", "No"))

mm.diff.child<- mm.diff.child%>% 
  mutate(statsig= if_else(p < 0.05, "Yes", "No"))

mm.diff.party<- mm.diff.party%>% 
  mutate(statsig= if_else(p < 0.05, "Yes", "No"))

SG1<-plot(mm.diff.gender, feature_headers = FALSE) +
  aes(color = statsig)+
  ggplot2::facet_wrap(~feature, ncol = 1L, labeller=labeller(feature = labels),
                      scales = "free_y", strip.position = "right") +
  theme_ggdist()+
  geom_label(aes(label=sprintf("%0.2f", round(estimate, digits = 2))), nudge_x = -.12, size=2.5)+
  xlim(-.3,.25)+
  scale_color_manual(values = c("black","red2"))+
  labs(subtitle="Gender differences")+
  xlab("Women - Men")+
  theme(
    axis.title.y =element_blank(),
    legend.position ="none",
    strip.text = element_blank(),
    legend.text = element_text(face="bold"),
    plot.subtitle = element_text(face="bold"),
    axis.title.x =element_text(face="bold"),
    axis.text.y = element_text(face="bold"),
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.title=element_blank()) +
  geom_vline(xintercept = 0.50, linetype="dashed", colour="gray33")

SG2<-plot(mm.diff.LGB, feature_headers = FALSE) +
  aes(color = statsig)+
  ggplot2::facet_wrap(~feature, ncol = 1L, labeller=labeller(feature = labels),
                      scales = "free_y", strip.position = "right") +
  theme_ggdist()+
  geom_label(aes(label=sprintf("%0.2f", round(estimate, digits = 2))), nudge_x = -.12, size=2.5)+
  xlim(-.3,.25)+
  scale_color_manual(values = c("black","red2"))+
  labs(subtitle="LGBT+ differences")+
  xlab("LGBT+ - Heterosexual")+
  theme(
    axis.title.y =element_blank(),
    legend.position ="none",
    strip.text = element_blank(),
    plot.subtitle = element_text(face="bold"),
    axis.title.x=element_text(face="bold"),
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.title=element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),) +
  geom_vline(xintercept = 0.50, linetype="dashed", colour="gray33")


SG3<-plot(mm.diff.child, feature_headers = FALSE) +
  aes(color = statsig)+
  ggplot2::facet_wrap(~feature, ncol = 1L, labeller=labeller(feature = labels),
                      scales = "free_y", strip.position = "right") +
  theme_ggdist()+
  geom_label(aes(label=sprintf("%0.2f", round(estimate, digits = 2))), nudge_x = -.12, size=2.5)+
  xlim(-.3,.25)+
  scale_color_manual(values = c("black","red2"))+
  labs(subtitle="Parental differences")+
  xlab("Parents - Not parents")+
  theme(
    axis.title.y =element_blank(),
    legend.position ="none",
    strip.text = element_text(face="bold"),
    legend.text = element_text(face="bold"),
    plot.subtitle = element_text(face="bold"),
    axis.title.x=element_text(face="bold"),
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.title=element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),) +
  geom_vline(xintercept = 0.50, linetype="dashed", colour="gray33")

SG4<-plot(mm.diff.party, feature_headers = FALSE) +
  aes(color = statsig)+
  ggplot2::facet_wrap(~feature, ncol = 1L, labeller=labeller(feature = labels),
                      scales = "free_y", strip.position = "right") +
  theme_ggdist()+
  geom_label(aes(label=sprintf("%0.2f", round(estimate, digits = 2))), nudge_x = -.12, size=2.5)+
  xlim(-.3,.25)+
  scale_color_manual(values = c("black","red2"))+
  labs(subtitle="Partisan differences")+
  xlab("Conservative - Labour")+
  theme(
    axis.title.y =element_blank(),
    legend.position ="none",
    strip.text = element_text(face="bold"),
    legend.text = element_text(face="bold"),
    plot.subtitle = element_text(face="bold"),
    axis.title.x=element_text(face="bold"),
    legend.background = element_blank(),
    legend.key = element_blank(),
    legend.title=element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),) +
  geom_vline(xintercept = 0.50, linetype="dashed", colour="gray33")


SG1+SG2+SG4
ggsave("figures/Figure5.eps", height=24, width=22, units="cm", dpi=900)

